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Abstract. A toad in the Rhinella granulosa group has been recognized as present on Trinidad since 1933. In 1965, the Trini¬ 
dadian population was described as a subspecies of Bufo granulosus, B. g. beehei. It has its type locality on the island and 
was eventually raised to species status as B. beehei (Beebe’s toad). Recently Beebe’s toad was synonymized with Rhinella 
humboldti, a species with a type locality in the Magdalena Valley of western Colombia. The Magdalena Valley is separated 
from the Orinoco Basin by the Eastern and Merida Cordilleras. These ranges have peaks that exceed 5,000 m and an al¬ 
most continuous altitude at about 3,000 m. Here we examine the morphology, advertisement calls, and mtDNA from sev¬ 
eral populations of these lowland toads to test whether the western Colombian R. humboldti and the Orinoco-Trinidad 
R. beebei are conspecific and form a single taxon that occurs on both sides of the Andes. The morphological, molecular, 
and advertisement call analyses suggest that R. humboldti and R. beebei are distinct taxa composed of independent evolv¬ 
ing lineages. Rhinella beebei is therefore resurrected from the synonymy of R. humboldti for the Trinidad and some of the 
adjacent mainland Orinoco populations in both Venezuela and Colombia. This increases the number of described species 
in the clade to fourteen. Rhinella humboldti and its sister R. centralis (Panama) are the only members of the R. granulosa 
group to occur west of the Andes, and our molecular results suggest the TMRCA for R. beebei and R. humboldti at about 
9 Mya, a time when the Eastern Cordillera was much lower in altitude than it is today and the Merida Cordillera was in 
its early stages of formation. 
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Introduction 

The timing of the Andean uplifts and related processes are 
critical to understanding Neotropical biogeography. Fos¬ 
sil fish of the genera Brachyplatystoma, Lepidosiren, Ara- 
paima, Phractocephalus, and Colossoma and fossils of the 
Colombian matamata turtle, Chelus colombiana, suggest 
that an ancient Amazonian-Orinoco fauna once occurred 
on both sides of the Andes (Lundberg 2005, Cadena et 
al. 2008). Prior to the existence of the Andean Eastern 
Cordillera, central Colombia was part of the palaeo-Ama- 
zon-Orinoco system (Hoorn et al. 1995, 2010, Lundberg 
et al. 1998). The north-flowing palaeo-Amazon-Orinoco 


drained the basin east of the developing Andean moun¬ 
tain chain from Bolivia to Venezuela. However, during or 
after the uplift of the Eastern and Northern Andes, the 
new drainage divide isolated the trans-Andean Magdale¬ 
na Valley, and much of the Amazonian and Orinoco fau¬ 
na disappeared from the region (Lundberg 2005). Today, 
the Eastern Cordillera (Colombia) and the Cordillera de 
Merida (Venezuela) form barriers between the Orinoco 
drainage to the east and the Magdalena and Maracaibo 
drainages to the west. These barriers have peaks that ex¬ 
ceed 5,700 m and altitudes that are almost continuous at 
3,000 m above sea level. Resolving estimates of barrier 
uplift and river drainage shifts with molecular-clock es- 
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timates of genetic divergence times for various Andean 
and Amazonian populations is important for understand¬ 
ing South Americas spectacular biodiversity (Javadi et al. 
2011). 

Two species of toads in the genus Rhinella are present 
in the Republic of Trinidad and Tobago. The marine toad, 
R. marina, is a large, widespread species known to reach at 
least 240 mm in body length and possibly the largest and 
best-studied toad in the world. It is locally known as the 
crapaud and has played a role in the islands’ ecosystems as 
well as the culture of Trinidad and Tobago (Kenny 1969, 
Murphy 1997). A second, smaller (< 60 mm) species of 
Rhinella is also present on Trinidad but not Tobago. This 
toad’s presence on the island has been known since at least 
1933 when Parker (1934) included it in his list of Trinida¬ 
dian anurans. It has been known to science as Bufo granu¬ 
losus, B. granulosus beebei, B. beebei, Rhinella beebei, and 
most recently as R. humboldti (Kenny 1969, Murphy 1997, 
Narvaes & Rodriquez 2009). 

A short historical summary of the Rhinella granulosa 
group nomenclature 

Spix (1824) described Bufo granulosus from Brazil, 
Dumeril & Bibron (1841) described B. dorbignyi from 
Uruguay, and Parker (1935) considered dorbignyi a sub¬ 
species of granulosa. Following the trends of consolidat¬ 
ing populations and names and recognizing polytypic spe¬ 
cies, Muller & Hellmich (1936) described B. g. major and 
considered B. granulosus to be composed of three subspe¬ 
cies: B. g. granulosus, B. g. dorbignyi, B. g. major. However, 
other authors suggested the subspecies might be distinct 
at species level (Schmidt & Inger 1951). Myers & Car¬ 
valho (1952) described B. pygmaeus from Rio de Janeiro, a 
species they believed to be closely related to B. granulosus. 
Gallardo (1957) described B. g. fernandezae from Buenos 
Aires, Argentina, and later Gallardo (1965) reviewed the 
group, recognizing 14 subspecies of B. granulosus that were 
distributed from the north coast of South America (includ¬ 
ing Trinidad and the Isla de Margarita) to Uruguay. Adding 
nine new taxa to those already described, he wrote, 

.. one can see in B. granulosus a clear distribution of the 
subspecies according to the hydrographical systems, a dis¬ 
tinct subspecies belonging to the basin of each of the follow¬ 
ing rivers: Magdalena, Orinoco, Upper and Middle Amazon, 
Tocantins, and Araguaia, Sao Francisco, Rio de la Plata...” 

Gallardo went on to comment that other distinct taxa 
were found in Guyana ( Rhinella merianae ), on the Isla de 
Margarita (R. barbouri ), and that the Trinidad population 
was the same as the one represented in the Orinoco drain¬ 
age (R. beebei). He designated a specimen from Trinidad as 
the holotype for Bufo granulosus beebei, and an additional 
13 specimens as paratypes, all with the same locality data 
(Churchill-Roosevelt Highway, Trinidad) 

Since Gallardo’s work, the polytypic species concept 
has been replaced with a lineage-based species concept 
(de Queiroz 1998) and advances in molecular technology 


have led to the resurrection of many old names and the dis¬ 
covery of many new species (Santos et al. 2009). The fam¬ 
ily Bufonidae was revised by Frost et al. (2006) and found 
to be a worldwide clade with 35 genera. The cosmopolitan, 
polyphyletic genus Bufo was subsequently divided into sev¬ 
eral genera, most of which follow the principles of bioge¬ 
ography. Rhinella Fitzinger, 1826 was the oldest available 
name for the clade of Neotropical toads that included both 
the R. granulosa and R. marina groups. 

Narvaes & Rodrigues (2009) re-examined the Rhinel¬ 
la granulosa group, using a Canonical Discriminant Anal¬ 
ysis and raised many of Gallardo’s subspecies of granu¬ 
losa to species level, synonymized others, and described a 
new species (R. centralis) from Panama, based on speci¬ 
mens Gallardo had considered to be Bufo granulosus 
humboldti. The end result was 12 species in the R. granulosa 
group. They viewed members of the granulosa group as not 
necessarily associated with drainages, but with open habi¬ 
tats with high rainfall, noting that most species (R. azarai, 
R. bergi, R. centralis, R. dorbignyi, R. fernandezae, R. hum¬ 
boldti, R. major, R. mirandaribeiroi, R. pygmaea) are found 
below 300 m above sea level. This is relevant to the Orino¬ 
co Basin-Trinidad populations because they regarded B. g. 
beebei Gallardo as a junior synonym of B. g. humboldti 
Gallardo. 

On the basis of biogeography alone, the synonymous sta¬ 
tus of R. beebei in R. humboldti is surprising. The straight- 
line distance between Girardot, Colombia and Trinidad is 
about 1,580 km. Keeping in mind that toads of the granulo¬ 
sa group are for the largest part lowland, savannah species, 
the topography between the type locality for humboldti and 
Trinidad includes the Cordillera Oriental, the Serrania de 
Perija, the Cordillera de Merida, the Coastal Cordillera of 
Venezuela, as well as a thousand kilometres or more of sa¬ 
vannah. Thus, the evolution of the South American land¬ 
scape has supplied multiple historic opportunities for cor¬ 
ridors as well as multiple barriers to gene flow before, dur¬ 
ing, and after the formation of the Andes. 

Gallardo (1965:114) described Bufo granulosus beebei 
based on AMNH 55774 (as well as 11 other specimens with 
the same locality data). The toads were collected along the 
Churchill-Roosevelt Highway in Trinidad at an altitude 
of about 30 m. Three pages later he described B. g. hum¬ 
boldti based upon a male specimen (MCZ 24882) from the 
municipalities of Girardot and Gualanday, Department 
of Tolima, Colombia (~4 °i6”N, 74°8i’W, altitude -326 m 
a.s.l.) and a female paratype MCZ 8978 from the munici¬ 
pality of Fundacion, Department of Magdalena, Colombia 
(~io°29’N, 74°i2’W, altitude ~6o m a.s.l.): both localities 
are in the Rio Magdalena drainage basin, but from oppo¬ 
site ends of the basin and more than 700 km apart. 

Authors writing about this toad in Trinidad prior to 
Gallardo’s work called it Bufo granulosus (Parker 1934). 
Post 1965, the commonly used name was Bufo granulosus 
beebei (Kenny 1969) until Rivero et al. (1986) elevated the 
subspecies to species level. Murphy (1997) followed using 
the name B. beebei and considered the toad an Orinoco/ 
Trinidad endemic. 
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Narvaes & Rodrigues (2009:43) chose the name hum- 
boldti over beebei “...because it is associated with continen¬ 
tal populations while beebei was named after insular speci¬ 
mens obtained at Trinidad and Tobago.” The ICZN code 
(Ride 1999) does not recognize page priority, thus this is 
an acceptable decision, assuming the two names actually 
represent a single operational taxonomic unit. 

Guerra et al. (2011) examined the advertisement calls 
of seven species within the granulosa group and found the 
calls consisted of long trills, composed of notes with a vari¬ 
able pulse number (2-8). Torres-Suarez & Vargas-Sali- 
nas (2013) reported on the call of R. humboldti from the 
Reserve of Wisirare, municipality of Orocue, Department 
of Casanare, Colombia (4°54’4o”N, 7i°26’i2”W, altitude 
126 m a.s.l.). This location lies within the Orinoco drain¬ 
age, and the toads call has four pulses per note. 

All GenBank sequences from Trinidad and the east 
side of the Andes have been labelled Rhinella hum¬ 
boldti (KP685211.1, KP685210.1, KP685174.1, KP685173.1, 
KP685153.1, KP685131.1, KP685099.1, KP685058.1, 

KP685025.1, KP685025.1, KP684965.1, KP684964.1, 

KP149216.1, KP149488.1, KP149473.1, DQ158276.1, 

DQ158434.1, DQ158358.1, EF532287.1, EF532269.1, 

EF532251.1). The molecular analysis and phylogeny of the 
R. granulosa group by Pereyra et al. (2015) found it con¬ 
sisted of 13 species distributed throughout open habitats 
of South America and Panama. They performed separate 
phylogenetic analyses under direct optimisation (DO) of 
nuclear and mitochondrial sequences and recovered the 
R. granulosa group as monophyletic. However, they found 
a topological incongruence that they explained as the re¬ 
sult of multiple processes of hybridisation and introgres- 
sion. All the sequences they used from R. humboldti were 
from the east side of the Andes. None were from the type 
locality of R. humboldti, nor from the west side of the An¬ 
des. 


Aims of this study 

Here we examine toads of the R. granulosa group from 
Trinidad, Venezuela, and Colombia to resolve the correct 
name for the island and Orinoco Basin populations and at¬ 
tempt to identify whether the species present on Trinidad is 
a widespread species found across northern South America 
on both sides of the Andes as hypothesized by Narvaes & 
Rodriquez (2009) or an Orinoco Basin/Trinidad endemic 
as suggested by Gallardo (1965) and Murphy (1997). 

Methods and materials 

Molecular methods 

DNA extraction, purification and amplification protocols 
follow Jowers et al. (2014). DNA was extracted from three 
liver samples of Rhinella humboldti adult males from Trini¬ 
dad (Republic of Trinidad and Tobago, close to Trincity 
Central Road, (UWIZM. 2012.27.72.1-3) and from three Co¬ 


lombian individuals from two localities east of the Colom¬ 
bian Andes (ICN 55784), Department of Casanare, Mu¬ 
nicipality of Trinidad, Vereda La Canada, Finca La Palmi- 
ta (5°i9’n,4”N, 7i°2o’5i”W) (ICN 55776), Department of 
Casanare, Municipality of Paz de Ariporo, Vereda La Co- 
lombina, Finca El porvenir (6°2’36.5”N, 7i°5’34.2”W) and 
west (CZUT-A 1717), Department of Tolima, Municipality 
of Prado, Vereda El Caiman, Hydrolectric dam “Hidropra- 
do”, Eco-Hotel Palo de Agua, 03°44’55.6”N, 74°5o’23.5”W, 
altitude 376 m a.s.l.). The primers for the 12S rDNA: 
12SA 5 ’ - A A ACT GGG ATTAG ATAC C C C ACTAT- 3 ’, 12SB 
5’-GAGGGTGACGGGCGGTGTGT-3’ (Kocher et al. 
1989) and 16SL RDNA: 16SL s’-GCCTGTTTATCAAAAA- 
CAT-3’, 16SH 5’-CCGGTCTGAACTCAGATCACGT- 3’ 
(Palumbi 1996) amplified an approximately 430 and 550 
base-pair (bp) fraction, respectively. Templates were se¬ 
quenced on both strands, and the complementary reads 
were used to resolve rare, ambiguous base-calls in Se¬ 
quencer v. 4.9. After removing PCR primers and incom¬ 
plete terminal sequences, 347 and 481 bp of the 12S and 16S 
rDNA gene fraction, respectively, 836 bp total, were avail¬ 
able for analyses and aligned with the larger Rhinella data 
set (-2,450 bp in length 12S and 16S rDNA) from Pereyra 
et al. (2015) to improve support of internal branches. 

Preliminary analyses were performed with the complete 
mtDNA data set of the R. granulosa group. However, be¬ 
cause the aim of our study was not to assess the overall 
phylogenetic relationships within the group, but rather to 
assess the relationship between R. humboldti and R. bee¬ 
bei, we only included two haplotypes per species for the fi¬ 
nal analyses. As suggested by Pereyra et al. (2015), R. ber- 
nardoi mtDNA sequences showed signs of introgression, 
likely through hybridisation, and were therefore excluded 
from the analyses. Additionally, BLAST searches were con¬ 
ducted in GenBank and matches with high genetic affinity 
(~ 98%) were downloaded and included in the alignment, 
including shorter R. humboldti sequences of either the 12S 
or 16S rDNA gene fragments, but not available for both. Se¬ 
quences were aligned in Seaview v. 4.2.11 (Gouy et al. 2010) 
under ClustalW2 (Larkin et al. 2007) default settings. Ge¬ 
netic p-distances and standard errors (% ± SE) were calcu¬ 
lated using MEGA v. 6 (Tamura et al. 2011). 

We used the software BEAST v. 1.5.4 (Drummond & 
Rambaut 2007) to estimate coalescence times between 
species (TMRCAs). We ran the analyses under a specia- 
tion (Yule process) prior to using a strict molecular clock 
model. Inspection of trace files indicated that a strict-clock 
model could not be rejected (the parameter ‘Coefficient of 
Variance’ in relaxed clock analyses abuts zero), and thus, we 
report results based on this model only. Dates of divergence 
estimated for 16S rDNA-corrected distances adopted for 
other anurans have been estimated to correspond to 0.39- 
0.40% divergence per million years (Martinez-Solano 
et al. 2004, Manzanilla et al. 2007) and are in accord¬ 
ance with rates assumed for the 16S and 12S rRNA genes 
in other amphibians such as newts and true salamanders 
(Caccone et al. 1997, Veith et al. 1998), ranging between 
0.4 and 0.7% per million years and 0.44-0.54% for the 12S 
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rDNA in Rana (Sumida et al. 2000). Thus, in the absence 
of a reliable prior on the substitution rate or a well-dated 
calibration point, we specified a normal prior on the clock 
rate with a mean value of 0.0022 substitutions/site per line¬ 
age per million years for the combined 12S and 16S rDNA 
gene fragments. This calculation was obtained by consider¬ 
ing the mutation rate of each gene fragment of 12S and 16S 
rDNA and the length of each in the alignment. We do stress 
however that this is just an estimated approximation of di¬ 
vergence to assess and date possible observed patterns to 
allopatric events in the study area and hence caution needs 
to be exercised in interpreting any molecular clock esti¬ 
mates. Analyses were run for 100,000,000 generations, and 
posterior distributions of parameter estimates were visually 
inspected in Tracer v. 1.5 (Rambaut & Drummond 2007). 
Ten percent of the resulting topologies and parameter val¬ 
ues were discarded as burn-ins in TreeAnnotator v. 1.8.2. 
The most appropriate substitution model for the Bayesian 
Inference (BI) analysis was determined by the Bayesian 
Information Criterion (BIC) in jModeltest v. 2 (Posada 
2008). The tree was constructed using the BI optimality cri¬ 
teria under the best-fitting substitution model (GTR+I+G). 
MrBayes was used with default priors, Markov chain set¬ 
tings, and random starting trees. Each run consisted of four 
chains of 100,000,000 generations and was sampled every 
1,000th generation. A plateau was reached after few gen¬ 
erations with 25% of the trees resulting from the analyses 
discarded as burn-ins. Phylogenetic relationships between 
haplotypes for each locus were estimated using a Maximum 
Likelihood (ML) approach, as implemented in the software 
RAxML v. 7.0.4 (Silvestro & Michalak 2010), using the 
default settings. The 50% bootstrap consensus tree was 
built in PAUP v. 4 (Swofford 2002). All analyses were per¬ 
formed through the CIPRES platform (Miller et al. 2010). 


Morphological methods and call analysis 

We examined museum material related to R. humboldti 
and R. beebei. External morphological data were collected 
for 149 museum specimens (Appendix 1). Measurements 
taken included SVL (snout-vent length), femur length, 
tibia length, horizontal orbit diameter, head width, head 
height (depth), foot length, and were taken to the nearest 
0.1 mm using dial callipers. Statistical analyses were con¬ 
ducted with Microsoft Excel QI Macros, and a PCA was 
performed using Data. 

Figure 1 illustrates the locations from which we have 
examined specimens as well as the type localities for R. g. 
beebei (Gallardo), R. g. barbouri (Gallardo), and R. g. 
humboldti (Gallardo). Hyack et al. (2001) detail the 
problems associated with obtaining repeatable measure¬ 
ments on anurans, particularly preserved specimens that 
have not been carefully prepared. We agree that obtaining 
measurements that can be duplicated is at best problem¬ 
atic. Thus we use the morphometries to test our position 
only in the broadest sense and ask the reader to look main¬ 
ly at the ranges reported. 

Advertisement call analyses were recorded by Morley 
Read (Trinidad). Trans-Andean R. humboldti were re¬ 
corded at Department of Tolima, municipality of Coello 
by Sigifredo Calvijo Garzon and loaned by the Herpe¬ 
tology, Eco-physiology and Ecology Research Group from 
the Universidad del Tolima. The call from Hato, Guarico 
(Colombia), was available on the Internet (https://www. 
youtube.com/watch?v=9_vzpEdtLqY, accessed 30 October 
2014k Call recordings were analysed in PRAAT v. 6.0.14. 

The largest samples of specimens examined were from 
Puerto Ayacucho, Venezuela on the Orinoco River, about 
580 km southwest of Trinidad, and Camarata, Bolivar, 



Figure 1. Type localities (coloured wheels) and localities from which members of the Rhinella granulosa group were examined. Yellow 
localities are trans-Andean, red localities are Orinoco-Trinidad. 
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Venezuela about 480 km south of Trinidad. We also ex¬ 
amined smaller numbers of toads from Maracay and Ara- 
gua, Venezuela. The Colombian specimens examined were 
mostly from the west side of the Andes, but also included 
specimens from Meta on the east side of the Andes. 

Results 

Molecular results 

The three Trinidadian samples differed by 1 bp differ¬ 
ence, a transversion (G-A) in the 16S rDNA gene frag¬ 
ments between frogs 1 (UWIZM.2012.27.72.1) and 2 
(UWIZM.2012.27.72.2) and frog 3 (UWIZM2012.27.72.3) 
(A-G). Both eastern Colombian samples were identical 
with the only exception of one indel insertion. GenBank 
BLAST searches of the sequences obtained from the Trini¬ 
dad and east Colombian frogs matched Rhinella humbold- 
ti for 99%. The 12S rDNA sequence from Trinidad frog 3 
matched R. humboldti (KP685211 and DQ158434) for 100% 
and for 99% with all R. humboldti with the 16S rDNA gene 
fragment. The other Trinidad individual was identical to 
R. humboldti from Trinidad (DQ158434) for both gene frag¬ 
ments. The eastern Colombian R. humboldti matched for 
100% both 12S and 16S rDNA gene fragments of R. hum¬ 


boldti (KP685211 and KP149366). Lack of R. humboldti 12S 
rDNA sequences from west of the Andes resulted in a 99% 
match to its sister taxon, R. centralis. The 16S rDNA frac¬ 
tion matched R. humboldi for 99% (1 bp difference) and 
R. centralis with 6 and 7 bp differences. 

The best-fitting model for the BI tree was the GTR+I+G 
(-lnL=6655.89206, BIC=i3686.875302). All analyses re¬ 
covered a well-resolved monophyletic clade (BPP: 1) with 
R. humboldti (i.e., R. beebei) from Trinidad, Venezuela and 
Colombia and sister clade to R. merianae (BI BPP: 1.00, 
ML: 78%, Beast BPP: 1.00). More weakly supported in the 
BI and RAxML analyses but recovered in all three analyses 
is the grouping of our sequenced frog from west Colom¬ 
bia with another R. humboldti from the same side of the 
Andes, ~ 500 bp of the 16S rDNA, which thus explains the 
moderate support, except in the Bayesian analyses with a 
0.96 BPP. This grouping forms a sister clade to R. centra¬ 
lis from Panama, with a 0.91 and 1.00 BPP in the BI and 
Beast tree respectively and with a 64% ML bootstrap sup¬ 
port. The R. humboldti (west Andes) + R. centralis sister 
relationship to R. humboldti (east Andes) + R. merianae is 
well supported in all analyses (BPP: 0.81, in the BI tree and 
1.00 BPP in the Bayesian analyses) (Fig. 2). 

Our results suggest an old divergence between R. hum¬ 
boldti from each side of the Andes (mean value for the TM- 
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Figure 2. Best ML tree for the 12S and 16S rDNA data set of the Rhinella granulosus species group. Taxonomic assignment follows R. 
beebei (east of the Andes) and R. humboldti (west of the Andes). Values on and under nodes are posterior probabilities recovered from 
Bayesian inference analyses and the 50% majority rule consensus tree, respectively. The three specimens with asterisks were sampled 
and sequenced for this study. Rhinella beebei and R. humboldti are shown in red. 
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RCA: 9 million years, with the 95% HPD interval ranging 
between 7.3 and 11 Mya). However, it is important to keep 
in mind that these estimates represent divergence times 
between haplotypes and thus necessarily precede popula¬ 
tion divergence times (although the magnitude of this off¬ 
set between haplotype and population divergence is un¬ 
known). The topologies resulting from BEAST analyses are 
fully concordant with those resulting from MrBayes and 
RAxML analyses, including support for the positioning of 
R. humboldti and R. beebei with their likely sister species, 
R. merianae and R. centralis (Figs 2 and 3). 

Genetic divergences between Trinidad and Venezuela 
(0.60% ± 0.02) were similar to the differences encountered 
between eastern Colombia (0.60% ± 0.07) and much high¬ 
er than within Trinidad localities (0.31% ± 0.14). Within 
R. beebei, the genetic divergence was relatively low (0.5% 
± 0.06). Comparison within R. humboldti was not possi¬ 
ble due to the lack of the 12S rDNA gene fragment (from 
KP149421, 16S rDNA: 0.21%). The combined 12S and 16S 
rDNA divergence between R. beebei and R. humboldti was 
1.44% ± 0.06, i.e., slightly higher than that between R. cen¬ 
tralis and R. humboldti (1.23% ± 0.26) and lower than be¬ 
tween R. centralis andR. beebei (2.39% ± 0.21) and between 
R. merianae and R. beebei (2.71% ± 0.18) Because of the 3% 
cut-off divergence estimated for taxonomical species de¬ 
limitations inferred from the 16S rDNA gene (Vieites et 
al. 2009), we here report on genetic divergence based only 
on this gene for comparative reasons: R. beebei andR. hum¬ 
boldti (1.98% ± 0.06); R. centralis and R. humboldti (1.0% ± 


0.09); R. centralis and R. beebei (2.4% ± 0.01); and R. meri¬ 
anae and R. beebei (2.7% ± 0.04). 

Morphological and call analysis results 

The examination of toads from Trinidad, several locations 
in Venezuela, and several locations in Colombia suggests 
multiple morphs to exist within R. humboldti. Two distinct 
morphs at Puerto Ayacucho suggest possibly sympatric 
species. 

Tibia/femur and femur/SVL ratios tested with two sam¬ 
ple t-tests assuming equal variance could not distinguish 
Trinidad and Venezuelan Camarata males (p > 0.05). How¬ 
ever, the tests readily distinguished the femur/SVL and tib¬ 
ia/femur ratios of the Puerto Ayacucho, the Trinidad-Ori- 
noco Basin populations, and the trans-Andean population 
(p < 0.05) from each other (Fig. 4). These traits and three 
others were tested with a PCA (details in Appendices 2-4). 

Advertisement calls: Figure 5 compares spectrograms 
and oscillograms for calls from a Trinidad specimen, an 
Orinoco Basin specimen (Venezuela), and a Trans-An¬ 
dean Colombian specimen. Frequency information was as 
follows: Trinidad - mean frequency 3,150 Hz, lowest fre¬ 
quency 2,790 Hz, maximum frequency 3,540; Orinoco Ba¬ 
sin - mean frequency 2,815 Hz, lowest frequency 2,400 Hz, 
maximum frequency 3,350 Hz; Colombia - mean frequen¬ 
cy 2,750 Hz, lowest frequency = 2,270 Hz, maximum fre¬ 
quency = 3,380 Hz. 



R. beebei * UWIZM.2012.27.72.3 Trinidad 
R. beebei (KP685210) Venezuela 

R. beebei * 1CN55776/8 (KP685211-KP149366) East Colombia 
R. merianae (KP685221) Brazil 
R. centralis (KP68196) Panama 
R. humboldti (KP149421) West Colombia 
R. humboldti * CZUT-A 1717 West Colombia 
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Figure 3. Bayesian inference tree obtained from BEAST of all major Rhinella granulosus group lineages. Naming follows R. beebei 
(east of the Andes) and R. humboldti (west of the Andes). Values at nodes are posterior probabilities recovered from the analyses; 
95% HPD intervals are shown as bars and ranges. Red values are the mean time since the most recent common ancestor (TMRCA, 
in million years). Species denoted with an asterisk were specimens sampled and sequenced for this study. The drawing in the tree 
represents the Andes mountain range. 
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Figure 4. A PCA of ratios calculated from members of the Rhinella humboldti group. Three distinct groups: green triangles and purple 
diamonds - Trinidad, Orinoco Basin, Bolivar Venezuela cluster; orange diamonds - Puerto Ayacucho, Venezuela cluster; and the blue 
squares and the red diamonds form the cis-Andean/trans-Andean Colombia cluster. Ratios included were femur:tibia, femur:SVL, 
orbit length:SVL, head width:head length, parotoid length:SVL, foot length:SVL The eigenvalues are given in Appendix 2. 


Calls from Trinidad and the Orinoco Basin (Venezue¬ 
la) specimens are similar in structure with four pulses per 
note, with the second and third the strongest, with similar 
amplitude and three notes per 1/10 of a second. The call 
from the Tolima population (Colombia) is divergent with 
only three pulses per note, with the second the strongest, 
and lower amplitude with four notes per 1/10 of a second. 
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Figure 5. Comparative spectrograms and oscillograms of the ad¬ 
vertisement call of Rhinella humboldti and Rhinella beebei from 
(A) Trinidad, (B) Guarico, Venezuela, and (C) Tolima, Colombia. 


Discussion 

The barriers to understanding speciation created by the 
polytypic species concept and subspecies was recognized 
and discussed by Cracraft (1983). Genetic divergences of 
the 16S rDNA gene between Rhinella beebei and R. hum¬ 
boldti (1.98%) are less than the 3.0% cut-off suggested by 
Vieites et al. (2009) to identify species level delimitation. 
However, those authors further qualify this value with a 
statement that the uncorrected pairwise genetic divergenc¬ 
es in the 16S rRNA gene to all other described species is 
in most cases 3%, but that it is sometimes only 1-2%. They 
furthermore state that the lower molecular differentiation 
can be augmented by a qualitative difference in advertise¬ 
ment calls, or a diagnostic difference in at least one mor¬ 
phological character known to be generally species-specif¬ 
ic. Moreover, applying the 3% cut-off threshold would re¬ 
sult in the lumping of R. beebei, R. humboldti, R. merianae 
and R. centralis in one species complex. We believe that the 
assignment of a certain accumulation of genetic divergence 
to establish taxonomic relationships should be treated with 
caution, as different taxa and evolutionary processes will 
affect the accumulation and or loss of genetic diversity dif¬ 
ferently. 

Our small sample size of vocalizations shows the cis-An- 
dean R. beebei has four pulses per note, while the trans-An- 
dean R. humboldti has three pulses per note. The number 
of pulses per note in the R. granulosus group was found to 
be species-specific by Guerra et al. (2011). 
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The morphological ratios used in our PCA, as well as 
such qualitative traits as the beaded ocular crests in R. bee- 
bei and the smooth ocular crests in R. humboldti, are ad¬ 
ditional morphological evidence (Fig. 6) for the independ¬ 
ent evolution of the populations on separate sides of the 
Andes and concur with the suggestions by Vieites et al. 
(2009). 

Given the molecular, morphological, and vocal evi¬ 
dence that the trans-Andean population is distinct from 
the cis-Andean populations, we remove Bufo granulosa 
beebei Gallardo from the synonymy of Rhinella hum¬ 
boldti (Gallardo), and recognize Rhinella beebei as the 
valid name for the Trinidad-Orinoco Basin member of the 
R. granulosa group. Our findings that two distinct lineages 
exist on either side of the Andes Mountains fully supports 
previous work (Guarnizo et al 2015) on the importance of 
the Andes mountain range as a key factor in allopatric spe- 
ciation in R. humboldti and other anurans. 

Furthermore, GenBank BLAST search results of west¬ 
ern Colombian R. humboldti matched the only sequenced 
R. humboldti from the west of the Andes, with a 99% af¬ 
finity and only 1 bp difference. Further evidence arises 
from the sister clade relationship to R. centralis from Pan¬ 
ama, the only other species within the R. granulosa group 
found west of the Andes. This finding was also support¬ 
ed by the lower genetic p-distances (12S and 16S rDNA) 



Figure 6. (A) Rhinella beebei from Trinidad. Photo by Joanna M. 
Smith; (B) R. humboldti from Antioquia, Colombia, municipality 
of Caucasia. Photo by Gustavo Gonzalez. 


from R. humboldti (1.23%) than between R. humboldti and 
R. beebei (1.44%). The sister relationship between R. meria- 
nae and R. beebei, both from the east side of the Andes, 
further corroborates the taxonomic distinction between 
R. humboldti and R. beebei. 

The timing of Andean uplifts and the establishment of 
a continuous topographic barrier in western South Amer¬ 
ica remain critical to understanding Neotropical biogeog¬ 
raphy. Horton (2014) noted that recent advances would 
allow independent geologic estimates of barrier uplift and 
river drainage shifts that can be compared with molecular- 
clock calculations of genetic divergence times for various 
Andean and Amazonian populations. 

Central to this discussion is Janzens (1967) hypothesis 
that mountain passes are “higher in the tropics.” Janzens 
climatic model predicts tropical mountain passes would be 
more effective barriers to a species’ dispersal than temper¬ 
ate-zone passes at equivalent altitudes. He predicted that 
tropical lowland organisms were more likely to encounter 
a mountain pass as a physical barrier to dispersal, which 
should in turn favour smaller distribution ranges and an 
increase in species turnover along altitudinal gradients 
(Janzen 1967, Ghalambor et al. 2006). 

Evidence suggesting the palaeoelevation of the Eastern 
Cordillera was lower in the past was discussed by Grego¬ 
ry-Wodzicki (2000). She reviewed the quantitative pal¬ 
aeoelevation estimates for the Central and Colombian An¬ 
des and suggested the Eastern Cordillera was 35-40% of 
the modern altitude (2,820 m above sea level) 5-4 Mya, 
which would correspond to about 1,000 m. Additionally, 
Javadi et al. (2011) notes that the Merida Andes were up¬ 
lifted in two phases. The first phase was in the Late Miocene 
(7-5 Mya) from the NW-SE convergence of the Panama 
Arc and western South America. The second phase in the 
Pliocene/Quaternary (5-1 Mya) resulted from the conver¬ 
gence of the Maracaibo Block and the Guyana Shield. Once 
these two ranges were established near or at their present 
altitudes, they formed an effective barrier that prevented 
lowland amphibians from moving across the mountains. 
Our approximate mean estimate of 9 Mya as the date of di¬ 
vergence is well before the 5 Mya date for the Eastern Cor¬ 
dillera uplift that would place the altitude at about 846 m, 
and is near the bottom of the time range for the uplift of the 
Meridian Andes (7-5 Mya). Divergence between R. beebei 
+ R. merianae and R. humboldti + R. centralis dates to the 
recession of the floodbasin system after the Miocene. The 
uplift of the Andes in the latter Miocene would have sepa¬ 
rated both clades. Similarly, a dendrobatid phylogeography 
(Santos et al. 2009) concluded that most extant Amazo¬ 
nian species dispersed into the Amazonian Basin after the 
Miocene floodbasin receded and immigration occurred 
on either side of the Andes around the Miocene-Pliocene 
transition. 

The combined barriers of the Eastern Cordillera and the 
Merida Andes with modern altitudes of 3,000-5,000 m are 
effective in preventing lowland species from moving east 
or west across these ranges. Rojas-Morales (2012) has 
provided the example of a forest-dwelling dipsadid snake, 
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Rhinobothryum, that apparently crossed the Andes from 
west to east and speciated in the Amazon basin to produce 
R. lentiginosum. In the case of the toads under discussion 
here, they likely moved west and speciated in the Magdale¬ 
na Valley ( R. humboldti) and adjacent Panama (R. centra¬ 
lis). 

The confusion of R. humboldti and R. beebei parallels 
the situation of four other anurans reported to occur in 
the humid lowlands on both sides of the Andes: Rhinel¬ 
la margaritifera, R. marina, Hypsiboas boans and Trachy- 
cephalus typhonius. There is genetic and morphological ev¬ 
idence that R. marina and Trachycephalus typhonius popu¬ 
lations on either side of the Andes represent separate spe¬ 
cies (Slade & Moritz 1998, Ron & Read 2011). Santos et 
al. (2015) assigned the populations of R. margaritifera from 
western Ecuador and Panama to R. alata and demonstrat¬ 
ed that the unusual distribution pattern of R. margaritifera 
on both sides of the Andes was an artefact of incorrectly 
defined species boundaries. 
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Appendix 1 

Specimens examined 

Colombia: Bolivar, Santa Rosa (io°27’N, 75°2o’W) ICN 2167, ICN 
2208; Boyaca, Puerto Boyaca Inspeccion, Puerto Romero, vere- 
da La Fiebre, finca La Barrilera (74°2o’29.47“N, 5 0 5o’i3.i2”W) 
ICN 38501, ICN 38504, ICN 38509; Vereda el Ocal, La Fiebre, 
Km. 29-30 carretera a Otanche (75°2o’39.73”N, 5°5o’25.32”W) 
ICN 38511; Puerto Boyaca, Inspeccion Puerto Romero, vereda 
a Fiebre, finca La Barrilera (74°2o’29.47”N, 5 0 5o’i3.i2”W) ICN 
44707; Quebrada La Cristalina, Km. 38-39 carretera a Otanche 
(74°i9T6.87”N, 5°49’40.7”W) ICN 45208-10; Campamento 
Technint. (74°29’33.9”N, 5 0 57’i7.32”W) ICN 45211; Boyaca Puerto 
Boyaca Vereda el Ocal, La Fiebre, Km. 29-30 carretera a Otanche 


(75°2o’39.73”N, 5°5o’25.32”W) ICN 45215-16; Caldas Norcasia via 
a la casa de maquinas, despues del Rio La Miel. (74°54’22.7”N, 
5°32’26.i5”W) ICN 40261; Caldas Samana km 15.6 carretera La 
Victoria-carretera central (74°56’40.8”N, 5°2i’36.27”W) ICN 
34711, 34713; Casanare, Paz de Ariporo Vd. La Colombina, Finca 
El porvenir (6°2’36.5”N, 7i°5’34.2”W) ICN 55776, ICN 55777, ICN 
55779-80; Trinidad,Vd. La Canada, Finca La Palmita (5 0 i9’n.4”N, 
7i°20.5i’W); ICN 55784; ICN 55781 (5 0 19 ’h. 4”N, 7i°2o’5i”W); 
ICN 55776 Paz de Ariporo Vd. La Colombina, Finca El Porvenir 
(6°2’36.5”N, 7i°5’34.2”W); ICN 55777-80. Cesar La Gloria Ciena- 
ga de Morales, poblacion Las Puntas (73°45’2”N, 8°32’42”W) 
ICN 37312; Huila, Campoalegre 75°i9’32”N, 2°4i’W) ICN 9367; 
Villavieja (75°i2’56.76”N, 3°i2’53.n”W) ICN 11819, 11821; Meta, 
San Martin Vereda La Castaneda, Cultivo de Palma Palma Sol 
S.A (3°3i’45”N, 73°32’3o”W) ICN 55782-3; ICN 55785 Cumaral, 
Hacienda La Cabana (4 °i 8’N, 73°2i’W); ICN 55786 Villavicen- 
cio. Vd. Santa Maria baja (4°i3’N, 73°38’W);Tolima Honda Ha¬ 
cienda Tupinamba (74°47’30.75”N, 5°8T5.i6”W) ICN 43838; ICN 
43948; Venadillo Finca Paloballo, km. 40 via Alvarado-Venadil- 
lo (74°56T7.8”N, 4°36’50 .i”W) ICN 43168-69; Venadillo 4 km al 
S de Venadillo (74°56’6.77”N, 4°4i’40.4”W) ICN 52093, 52098- 
99; Tolima (75°9T6.34”N, 4°5’33.6”W) ICN 17636; Mariquita, 
(~5°i4’N, 74°55’W) FMNH 81833-34, Tolima no specific locality 
FMNH 54182, ICN 3848, 17636, 43168-69, 43938, 52093, 52098- 
99; Trinidad: BMNH 22, 116, 546-47, 1547-48; FMNH 218784, 
251214. 

Venezuela: Aragua (~io° 07’N, 67°58’W) FMNH 69780, 
176430; Estado Bolivar, Camarata, (~5°38”N, 62 °i 8”W) BMNH 
149-151, 154-155, 157-159; Maracay (io° 14’N, 67°35’W) FMNH 
125405, 174359; FMNH Puerto Ayacucho (~27°58’N, 67°35’W) 
FMNH 175307-088, 175311, 175313-14, 175319, 175322, 175325, 

175327 . 175338 , 175342 . 


Appendix 2 

Results of the two sample t-tests assuming equal variances 

(A) Comparison of tibia/femur ratios from four populations of 
Rhinella beebeilhumboldti, alpha set at 0.05; (B) Comparison of 
femur/SVL ratios from four populations of Rhinella beebeil hum¬ 
boldti, alpha set at 0.05. Bold typeface indicates null hypothesis 
rejected. 


A 

Llanos, 

Venezuela 

Bolivar, 

Venezuela 

Tolima, 

Colombia 

Trinidad 

p=0.00, df=13 

p= 0 . 098 , df =14 

p=0.00, df=13 

llanos 


p=0.00, df =16 

p=0.00, df =15 

Bolivar 



p=0.00, df =15 

B 

Llanos, 

Venezuela 

Bolivar, 

Venezuela 

Tolima, 

Colombia 

Trinidad 

p=0.00, df=14 

p= 0 . 67 , df =14 

p= 0 . 114 , df =13 

llanos 


p=0.00, df =16 

p=0.00, df =14 

Bolivar 



p=0.00, df =15 
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Appendix 3 

The communalities and eigenvalues from the PCA 

Ratios derived of measurements included femur:tibia (f/t), femur:svl (f/SVL), orbit length:SVL ol/SVL, head width:head length (hw/ 
hi), parotoid length:SVL (pl/SVL), foot length:SVL (ft/SVL). 



ratios 

PC -1 

PC -2 

PC -3 

PC -4 

PC -5 

PC-6 

1 

f/t 

0.5190 

0.1867 

0.07682 

0.2052 

0.5145 

0.6189 

2 

f/SVL 

0.5051 

0.1328 

0.3075 

0.3123 

0.1364 

- 0.7188 

3 

ol/SVL 

- 0.4350 

0.1240 

0.09886 

0.8717 

- 0.1099 

0.1174 

4 

hw/hh 

0.1301 

0.8580 

- 0.09063 

- 0.1157 

- 0.4693 

0.07183 

5 

pl/SVL 

- 0.2675 

0.1199 

0.9099 

- 0.2632 

0.05241 

0.1190 

6 

ft/SVL 

0.4448 

- 0.4262 

0.2317 

0.1338 

- 0.6940 

0.2594 


eigenvalues 

3.160 

1.208 

0.8709 

0.4893 

0.1969 

0.07556 


Appendix 4 

Percentages of variation for each PC 
Numbers correspond to those in Appendix 3 . 


PC 

% variation 

total variation 

1 

53.42 

53.42 

2 

20.3 

73.73 

3 

14 

87.73 

4 

8.22 

95.95 

5 

3.31 

99.26 

6 

0.74 

100 
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